/** 
 *  @file
 *  @brief   小量drill 膜单元
 *  @author  wanghuan
 *  @date    2024.03.13
 *  @version  v1.0
 *  @par Copyright(c):  ShenGong  Simulation corporation
 */

#pragma once
#include "Utility/Algebra/Matrix.h"
#include "DataStructure/Input/Element/CQUAD4Data.h"
#include "DataStructure/Input/Property/ShellSecMaterialConstant.h"
#include "SGFem/Integration/IntegrationFormula/Gauss2D2X2P.h"
#include "SGFem/Integration/IntegrationFormula/Gauss2D1X1P.h"
#include "SGFem/Integration/ShapeFunction/Plane4NodeShpFunc.h"
#include "SGFem/Integration/ShapeFunction/Plane4NodeIncomShpFunc.h"
#include "MemElement.h"


/** @brief 小量Drill四边形膜单元
 * 
 * 通过增加一个小量提供膜的drill刚度
 */
class PublicElement_DECLSPEC  Quad4SmallDrillMem : public MemElement
{
public:
    Quad4SmallDrillMem () = default;
    /// @brief 计算和节点坐标无关的变量
    void Initialize ();

    /** @brief 计算应变矩阵B
     * 
     * @param [in] iGauss 当 isCenter = false, iGauss 表示积分点索引
     * @param [in] coord 单元节点坐标
     * @param [out] B    应变矩阵
     * @param [out] j 积分点雅克比的行列式值
     * @param [in]  isCenter 当希望计算单元中心应变矩阵，输入true
     */ 
    void ComputeStrainMatrix(std::size_t iGauss, const SG::Algebra::Matrix& coord, _OUT SG::Algebra::Matrix& B, _OUT SG::DataStructure::Common::Real& j, bool isCenter = false) const;

    /** @brief 计算非协调应变矩阵G
     * 
     * @param iGauss  gauss点索引
     * @param  j      gauss点Jacobian的行列式值
     * @param JT0     单元中心点Jacobian的转置
     * @param [out] G 非协调应变矩阵G
     */
    void ComputeIncomStrainMatrix (std::size_t iGauss, SG::DataStructure::Common::Real j, const SG::Algebra::Matrix& JT0, _OUT SG::Algebra::Matrix& G) const;

    /** @brief 计算单元积分点应变
     * 
     * @param[in] localDisp 单元局部坐标系下节点位移
     * @param[in] localCoord 单元局部坐标系下节点坐标
     * @param[in] Cm0       材料弹性矩阵
     * @param[in] thickness 积分点厚度列表
     * @param[in][out] integStrain 修改膜部分应变
     */
    void ComputeIntegStrain (
        const SG::Algebra::Matrix& localDisp, 
        const SG::Algebra::Matrix& localCoord, 
        const SG::Algebra::Matrix& Cm0,
        const SG::DataStructure::Common::Real* const thickness,
        SG::Algebra::Matrix& _OUT integStrain);

    /** @brief 计算单元中心点应变
     * 
     * @param[in] localDisp 单元局部坐标系下节点位移
     * @param[in] localCoord 单元局部坐标系下节点坐标
     * @param[in] thickness 中心点厚度
     * @param[in][out] integStrain 修改膜部分应变
     */
    void ComputeIntegStrain0 (
        const SG::Algebra::Matrix& localDisp, 
        const SG::Algebra::Matrix& localCoord,
        SG::Algebra::Matrix& _OUT integStrain);

    /** @brief 计算单元刚度矩阵
     * 
     * @param [in] coord    提供局部坐标系
     * @param [in] thickness 节点厚度数组
     * @param  SecConstant  截面属性，需要用到膜部分材料刚度 Cm（没有乘厚度 t ）
     * @param  T1           膜部分材料主轴->单元局部坐标系x轴转换矩阵
     * @param  transFalg    材料主轴与单元坐标系x轴平行标识，如果平行，则不需要材料坐标转换
     */ 
    SG::Algebra::Matrix ComputeStiffness 
    (
    const SG::Algebra::Matrix& coord, 
    const SG::DataStructure::Common::Real* const thickness,
    const SG::DataStructure::FEM::ShellSecMaterialConstant& SecConstant,
    const SG::Algebra::Matrix& T1,
    const bool & transFlag
    ) const;
   
private:
    SG::FEMSolver::Element::Gauss2D2X2P       m_integForm; ///< 积分格式
    SG::FEMSolver::Element::Plane4NodeShpFunc m_shpFunc; ///< 积分点上的形函数
    SG::FEMSolver::Element::Plane4NodeShpFunc m_shpFunc0; ///< 中心点形函数
    SG::FEMSolver::Element::Plane4NodeIncomShpFunc m_incomShpFunc; ///<非协调项
    SG::FEMSolver::Element::Gauss2D1X1P       m_integForm0; ///< 中心点积分格式

    static std::size_t s_nNode;
};